In [26]:
%pylab inline
import pydisdrometer
from matplotlib.colors import LogNorm
In [12]:
filename = '/Users/hard505/Dropbox/Projects/NetworkRetrieval/data/Disdrometer/IFloodS_APU01_2013_0502_dsd.txt' #Parsivel 05, March 2nd
In [13]:
dsd = pydisdrometer.read_parsivel_nasa_gv(filename, campaign='ifloods')
So at this point we have the drop size distribution read in. NASA strips out rainfall information though. Let's do some T-Matrix scattering though. This should take a little bit(up to a minute or so on my laptop depending on the file).
In [14]:
dsd.calculate_radar_parameters()
This assumes BC shape relationship, X band, 10C. You can pass in a new wavelength to change that. Let's plot some of the parameters, and then try to do something with the data.
In [15]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.plot(dsd.time/60.0, dsd.fields['Zh']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Reflectivity(dBZ)')
plt.xlim(5,24)
plt.subplot(2,2,2)
plt.plot(dsd.time/60.0, dsd.fields['Zdr']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Differential Reflectivity(dB)')
plt.xlim(5,24)
plt.subplot(2,2,3)
plt.plot(dsd.time/60.0, dsd.fields['Kdp']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Specific Differential Phase(deg/km)')
plt.xlim(5,24)
plt.subplot(2,2,4)
plt.pcolor(dsd.time/60.0, dsd.diameter, np.log10(dsd.Nd.T), vmin=0, vmax=6)
plt.xlabel('Time(hrs)')
plt.ylabel('Diameter')
plt.colorbar()
plt.ylim(0,5) #Zoom in some
plt.xlim(5,24)
plt.tight_layout()
plt.show()
Next let's estimate some microphysical parameterizations.
In [16]:
dsd.calculate_dsd_parameterization()
In [22]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.plot(dsd.time/60.0, dsd.fields['D0']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('$D_0$')
plt.xlim(5,24)
plt.subplot(2,2,2)
plt.plot(dsd.time/60.0, np.log10(dsd.fields['Nw']['data']))
plt.xlabel('Time(hrs)')
plt.ylabel('$log_{10}(N_w)$')
plt.xlim(5,24)
plt.subplot(2,2,3)
plt.plot(dsd.time/60.0, dsd.fields['Dmax']['data'])
plt.xlabel('Time(hrs)')
plt.ylabel('Shape parameters')
plt.xlim(5,24)
plt.subplot(2,2,4)
plt.pcolor(dsd.time/60.0, dsd.diameter, np.log10(dsd.Nd.T), vmin=0, vmax=6)
plt.xlabel('Time(hrs)')
plt.ylabel('Dmax')
plt.colorbar()
plt.ylim(0,5) #Zoom in some
plt.xlim(5,24)
plt.tight_layout()
plt.show()
At hour twenty we see some drops that are probably instrument inaccuracies. Let's look at the overall behavior of the parameters with respect to each other.
In [68]:
plt.figure(figsize=(12,12))
plt.subplot(2,2,1)
plt.hist2d(dsd.fields['D0']['data'], np.log10(dsd.fields['Nw']['data']),
range=((0,6),(0,6)), vmin=1, bins=50, norm=LogNorm())
plt.ylabel('log10(Nw)')
plt.xlabel('$D_0$')
plt.subplot(2,2,2)
plt.hist2d(dsd.fields['Zh']['data'], dsd.fields['Zdr']['data'],
range=((-10,40),(-1,1)), vmin=1, bins=50, norm=LogNorm())
plt.xlabel('$Z_h$')
plt.ylabel('$Z_{dr}$')
plt.subplot(2,2,3)
plt.hist2d(dsd.fields['Dm']['data'], dsd.fields['Zdr']['data'],
range=((0,2),(-1,1)), vmin=1, bins=50, norm=LogNorm())
plt.xlabel('$D_m$')
plt.ylabel('$Z_{dr}$')
plt.subplot(2,2,4)
plt.hist2d(dsd.fields['Dmax']['data'], dsd.fields['Zh']['data'],
range=((0,8),(-10,40)), vmin=1, bins=50, norm=LogNorm())
plt.ylabel('$Z_h$')
plt.xlabel('$Dmax$')
plt.tight_layout()
plt.show()
Let's take what might be another common use case. Finding areas with certain size drops. In this case everything above 4mm. There are a few ways to do this but let's just create an indicator function for now.
In [19]:
plot(dsd.time/60.0 , 0<sum(dsd.Nd[:,dsd.diameter>4], axis=1))
plt.xlim(5,24)
Out[19]:
Or maybe we just want a listing of the time steps where there are drops greater than 4mm
In [20]:
time_hrs = dsd.time/60.0
print(time_hrs[0<sum(dsd.Nd[:,dsd.diameter>4], axis=1)])
In [ ]: